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equations in spherical polar coordinates without any symmetry assumptions. Using a PIRK method 
we obtain stable simulations in three spatial dimensions without the need to regularize the origin or 
the axis. We perform and discuss a number of tests to assess the stability, accuracy and convergence 
of the code, namely weak gravitational waves, "hydro-without-hydro" evolutions of spherical and 
rotating relativistic stars in equilibrium, and single black holes. 
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I. INTRODUCTION 



The first announcements of successful binary black hole 
simulations [IH2] marked an important break-through 
in numerical relativity and triggered a burst of ac- 
tivity in the field. While most current simulations 
adopt some variation of the BSSN formulation to- 
gether with what have become "standard coordinates" 
(namely 1+log slicing [7] and the "Gamma-driver" con- 
dition 8 ), different implementations differ in many de- 
tails. Most current, three-dimensional numerical relativ- 
ity codes share one feature, though, namely Cartesian co- 
ordinates. While Cartesian coordinates have many desir- 
able properties, there are applications, for example grav- 
itational collapse and supernova calculations, for which 
spherical polar coordinates would be better suited. 1 

Implementing a numerical relativity code in spherical 
polar coordinates poses several challenges. The first chal- 
lenge lies in the equations themselves. The original ver- 
sion of the BSSN formulation, for example, explicitly as- 
sumes Cartesian coordinates (by assuming that the deter- 
minant of the conformally related metric be one). This 
issue has been resolved by Brown |10) . who introduced 
a covariant formulation of the BSSN equations that is 
well-suited for curvilinear coordinate systems (compare 

Em)- 

Another challenge is introduced by the coordinate sin- 



Cartesian or spherical polar coordinates are not the only two 
possibilities, of course. In particular, multi-patch applications, 
combining some of the advantages of both, may present an at- 
tractive alternative at least for some applications (see, e.g., [9] 
for an implementation in numerical relativity). 



gularities at the origin and the axis, which introduce sin- 
gular terms into the equations. Although the regularity 
of the metric ensures that, analytically, these terms can- 
cel exactly, this is not necessarily the case in numerical 
applications, and special care has to be taken in order to 
avoid numerical instabilities. 

Several methods have been proposed to enforce reg- 
ularity in curvilinear coordinates. One possible ap- 
proach is to rely on a specific gauge, e.g. polar-areal 
gauge [HI [T3], together with a suitable choice of the 
dynamical variables. Numerous different such methods 
have been implemented in spherical symmetry (see, e.g., 
[14] for an overview); examples in axisymmetry include 
[THEEH1I. This approach has some clear limitations. It 
is not obvious how to generalize these methods to relax 
the assumption of axisymmetry; moreover the restriction 
of the gauge freedom prevents adoption of the "standard 
gauge" that proved to be successful in evolutions with 
the BSSN formulation. 

An alternative method is to apply a regularization pro- 
cedure, by which both the appropriate parity regular- 
ity conditions and local flatness are enforced in order to 
achieve the desired regularity of the evolution equations 
(see [TTl - f2"5] for examples). Typically, these schemes in- 
volve the introduction of auxiliary variables as well as 
finding evolution equations for these variables. The re- 
sulting schemes are quite cumbersome, which may ex- 
plain why, to the best of our knowledge, no such scheme 
has been implemented without any symmetry assump- 
tions. 

In yet an alternative approach, Cordero-Carrion et.al. 
[2~3] recently adopted a partially implicit Runge-Kutta 
(PIRK) method (see also [25]) to evolve the hyperbolic, 
wave-like equations in the Fully Constrained formula- 
tion of the Einstein equations (see [H]). Essentially, 
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PIRK methods evolve regular terms in the evolution 
equations explicitly, and then use these updated values 
to ev olve singular terms implicitly (see [55] and Section 
III A below for details). Following this success, Montero 
& Cordero-Carrion [27] . assuming spherical symmetry, 
applied a second-order PIRK method to the full set of 
the BSSN Einstein equations in curvilinear coordinates, 
and produced the first successful numerical simulations 
of vacuum and non- vacuum spacetimes using the covari- 
ant BSSN formulation in spherical coordinates without 
the need for a regularization algorithm at the origin (or 
without performing a spherical reduction of the equa- 
tions, compare [2"51 |2"§] V 

In this paper we present a new numerical code that 
solves the BSSN equations in three-dimensional spherical 
polar coordinates without any symmetry assumptions. 
The code uses a second-order PIRK method to integrate 
the evolution equations in time. This approach has the 
additional advantage that it imposes no restriction on 
the gauge choice. We consider a number of test cases 
to demonstrate that it is possible to obtain stable and 
robust evolutions of axisymmetric and non-axisymmetric 
spacetimes without any special treatment at the origin or 
the axis. 

The paper is organized as follows. In Section |ll] we 
present the basic equations; we will review the covari- 
ant formulation of the BSSN equations, and will then 
specialize to spherical polar coordinates. In Section III 



we will briefly review PIRK methods and will then de- 
scribe other specifics of our numerical implementation. 



In Section IV we present numerical examples, namely 
weak gravitational waves, "hydro-without-hydro" simu- 
lations of static and rotating relativistic stars, and single 
black holes. Finally we summarize and discuss our find- 
ings in Section[Vj We also include two appendices; in Ap- 
pendix[A]we describe an analytical form of the flat metric 
in spherical polar coordinates that provides a useful test 
of the numerical implementation of curvature quantities, 
while in Appendix [B] we list the specific source terms for 
our PIRK method applied to the BSSN equations. 

Throughout this paper we use geometrized units in 
which c = G = 1. Indices a, b, . . . denote spacetime in- 
dices, while i, j, ■ ■ ■ represent spatial indices. 



unity, which completely determines the conformal factor. 
This approach is suitable when Cartesian coordinates are 
used, but not in more general coordinate systems. We 
will pose a different condition on 7 below, but note al- 
ready that 
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(7/7) 



1/3 



(2) 



The advantage of this approach is that all quantities in 
this formalism may be treated as tensors of weight zero 
(see also [H]). We also denote 



K, 



1 



(3) 



as the conformally rescaled extrinsic curvature. Slightly 
departing from Brown's approach we assume this quan- 
tity to be trace-free, while Brown allows A^j to have a 
non-zero trace. In the above expression is the phys- 
ical extrinsic curvature and K — j t3 Kij its trace. 

Introducing a background connection Tj k (compare 
[TT] ) we now define 



r 



(4) 



which, unlike the two connections themselves, transform 
as a tensor field. We also define the trace of these vari- 
ables as 



Ar l = f' fc Ar!- 



■jk- 



(•5) 



It is not necessary for the background connection to be 



associated with any metric. In Section II B below we will 
specialize to applications in spherical polar coordinates 
and hence will assume that the T l j k are associated with 
the flat metric in spherical polar coordinates. This as- 
sumption affects the equations in the remainder of this 
Section in only one way, namely, we will assume that the 
Riemann tensor associated with the connection T l - k van- 
ishes, as is appropriate when the background metric is 
flat. 

Finally, we define the connection vector A* as a new 
set of independent variables that are equal to the Ar l 
when the constraint 



II. BASIC EQUATIONS 

A. The BSSN equations in covariant form 

We adopt Brown's covariant form [TO] of the BSSN 
formulation @HB]. I n particular, we write the conformally 
related spatial metric 7y as 

7ij e T»j 1 (1) 

where 7^ is the physical spatial metric, and a con- 
formal factor. In the original BSSN formulation the de- 
terminant 7 of the conformally related metric is fixed to 



C ee A 1 - AF = (6) 

holds. The vector A* plays the role of the "conformal con- 
nection functions" F in the original BSSN formulation, 
but, unlike the f l , the A* transform as a rank-1 tensor 
of weight zero (compare exercise 11.3 in [2]). In the 
following we will evolve the variables A 1 as independent 
variables, satisfying their own evolution equation. 

In order to determine the conformal factor we spec- 
ify the time evolution of the determinant of the conformal 
metric. In this paper we adopt Brown's "Lagrangian" 
choice 

da = 0. (7) 
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Defining 



d i =d t -C (- 



(8) 



where Cp denotes the Lie derivative along the shift vector 
f3 z , we then obtain the following set of evolution equations 



d±Aij 



%D k p k - 2aA. l 



(9a) 



-A l0 D k p k ~ 2aA lk A k j + aA i:j K 



+e 



-40 



- 2aD i D j i 
a{Ri - 8nSij] 



DiDja 



TF 



d ± K 



-4<p 



(D 2 a + 2D l aD l 



l -D k p - \aK 
6 6 

|x 2 + aAijRl 
+4?ra(p + S) 

f%V k l3 l + |aTVD^ + \&D^ 

-2A jk {5 l 3 d k a - Qa5 l 3 d k <t> - aAY) k ) 
4 

_ 3' 



-arf 3 djK — Myrrarf 13 Sj. 



(9b) 

(9c) 

,4>) 
(9d) 

(9e) 



(compare equations (21) in [TOj)- In the above equations, 
a is the lapse function, £>j denotes a covariant derivative 
that is built from the background connection T l - k (and 
hence, in our implementation, associated with the flat 
metric jij in spherical polar coordinates) and the super- 
script TF denotes the trace- free part. The matter sources 
p, Si, Sij and S denote the density, momentum density, 
stress, and the trace of the stress as observed by a normal 
observer, respectively, and are defined by 



p 


= n a n b T ab , 


(10a) 


Si 


= -J la n b T ab , 


(10b) 


Sij 




(10c) 


s 


— 7 


(lOd) 


n a 


= (-a, 0,0,0) 


(11) 



Here 



is the normal one-form on a spatial slice, and T ab is the 
stress-energy tensor. 

We compute the Ricci tensor Rij associated with 7^- 
from 



-~7*'© fc £,7<i + ^ k{ {D 3) K k + Ar k AT {lj)k 



"7 



; < (2Ait (i Ar i)ml + Ar^Ar 



mjt 



(12) 



In all of the above expressions we have omitted terms 
that include the Riemann tensor Rij k l associated with 



the connection r* fc , since these terms vanish for our case 
of a flat background. 

The Hamiltonian constraint takes the form 



H = -K 2 - AijAV +e- i4 '{R-SD i (j)D i (j>-d,D 2 



A6irp 



= 0, 



(13) 



while the momentum constraints can be written as 

M l = e-^(^-V ] {^A l3 )+(SA l3 d. J cp 
v V7 

9 .. . \ 

8irS l 



■jk 



= 0. 



(14) 



(see equations (16) and (17) in [TU]). 

We note that when 7 = 1 and T l - k = 0, which is suit- 
able for Cartesian coordinates, the above equations re- 
duce to the traditional BSSN equations. In the following, 
however, we will evaluate these equations in spherical po- 
lar coordinates. 

Before the above equations can be integrated, we have 
to specify coordinate conditions for the lapse a and the 
shift f3 l . Unless noted otherwise we will adopt a "non- 
advective" version of what has become the "standard 
gauge" in numerical relativity. Specifically, we use the 
"1+log" condition for the lapse [7] in the form 



d t a = —2aK, 



(15) 



and the "Gamma-driver" condition for the shift [5] in the 
form 



dtp 1 



B' 



d t B l = -d t V 
1 4 1 



(16a) 
(16b) 



(compare 27J). These (or similar) conditions play a 
key role in the "moving-puncture" approach to handling 
black hole singularities in numerical simulations. 



B. Implementation in spherical polar coordinates 

We now focus on spherical polar coordinates, and will 
assume that the Tj k are associated with the flat metric 
in spherical polar coordinates r, 0, and 



Hi = Va 







r 2 sin 2 1 



(17) 



Accordingly, the only non-vanishing components of the 
background connection are 



n e = -r 

tL. = -sin6»cos( 



f^ = -rsin 2 fl 

r% = r-i (18) 

ft, = cot0. 
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When implementing the above equations in spherical 
polar coordinates, care has to be taken that coordinate 
singularities do not spoil the numerical simulation. These 
singularities appear both at the origin, where r = 0, and 
on the axis where sin 9 = 0. Even for a simple scalar 
wave, appearances of inverse factors of r and sm9 in the 
Laplace operator can pose a challenge for a numerical 
implementation. In Section |TTl] below we discuss a PIRK 
method (see also [HI [57] ) that handles these singularities 
very effectively. 

An additional challenge in general relativity is that 
these inverse factors of r and sin# appear through the 
dynamical variables themselves. Components of the spa- 
tial metric, for example, scale with powers of r and sin 9, 
the inverse metric then scales with inverse powers of these 
quantities, and numerical error affecting these terms may 
easily spoil the numerical evolution. It is therefore im- 
portant to treat these appearances of r and sin 9 analyt- 
ically. We therefore factor out suitable powers of r and 
s'md from components of all tensorial objects. 2 

We start by writing the conformally related metric "fij 
as the sum of the flat background metric %j and a cor- 
rection eij (which is not assumed to be small), 



in terms of the coefficients hij. Direct calculation using 



(19) 



The flat metric 7^ is given by eq. (17), and we write the 
correction e ?1 in the form 



rh r g 
r 2 h e 
rsiu.9h rt p r 2 smOhg 



rh r g 



r sin 9h ri f, 
r 2 sin 9hg^ 
r 2 sin 2 6hd„t 



(20) 



We similarly rescale the extrinsic curvature Aij as 

ra r g r sin 9a r 



A, 




r 2 a ee 



r 2 sin 9ag ( f 1 
r 2 sin 2 9a 6d 



and the connection vector A 1 as 

X r 

A 1 = [ \ 9 /r 

A^/(rsin6») 



(21) 



(22) 



We treat the shift and B % similarly, and finally rewrite 
the evolution equations ^ for the coefficients hij , ay and 
X 1 etc. 

We can compute the connection coefficients Q from 

Ar;. fc = l -f l (v^u + T> k ^ - f>a Jk ) . (23) 

Since T>ijjk = we can compute the derivatives of the 
spatial metric 



(24) 



2 In an alternative approach, one could represent the metric in an 
orthonormal frame, so that the correct powers of r and sin 9 are 
absorbed in the unit vectors. 



the flat connection ( 18 1 yields 



D r r y rr 








= rh r g y r 






= rsm9h r< f, t r 






= r z hgg jr 






= r 2 sin 9hg c f > ^ r 






= r 2 sin 2 9h r f >r f > ^ r 






= hrr.e — 2h r g 




Vgjrg 


= r{h r g,g + h rr - hgg) 






= rsin9(h r4> . e - h g<j> ) 




£>e7ee 


= r 2 (hgg.g + 2h r g) 




£>070</> 


= r 2 s\Yi9(h S4> fi + hrf) 




'Del 4,4, 


= r 2 sin 2 Oh^g 






— h rr-( p — 2 sin 9h r( j) 




£>07re 


= r(h r g^ - cos 6h r<j> - sinfl/i^) 






= r sin Q{h T <l>,4> + si n 9h rr + cos 9h r g 


— sin 9h 


D^Jgg 


= r 2 (hgg^ - 2 COS 9hg<j ) ) 




2JV700 


= r 2 sin 0(/i0<4 <4 + sin 9h r g + cos 9hg 


9 — COS 9 


V4.I44 


= r 2 sin 2 0{h^^ + 2 sin 0h rt f, + 2 cos 6he<j>) 



The (flat) covariant derivative of the connection vector 
A 1 can similarly be expressed in terms of the A* as 

V r A r = d r X r 
V e A r = dg\ r -\ e 
V,\r = ^A r -sin6»A* 



V r A e 

VgA 6 



1 



d r \ e 



- (doX° + A") 



V^A 9 = - (<9 A e - cos 



(26) 



VgM 



r sm 
1 

r sin 
1 



sin9X r + cos9\ 6 



Using the above expressions, we can compute the Ricci 
tensor ( 12 ) as follows. In the first term on the right-hand 
side of ( 12 ) we write the second covariant derivative of 



jij as a sum of first partial derivatives of the quantities 
T>ijij and (flat) connection terms multiplying the PjTy, 



(27) 



We then insert the expressions ( 25 ) into the first term on 



the right-hand side and evaluate all derivatives explicitly, 
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so that these terms can be written in terms of second 
partial derivatives of the coefficients hij. Once this step 
has been completed, we add those remaining terms for 
which the flat background connection ( 18 ) is nonzero. 



The resulting equations are rather cumbersome, and 
it is easy to introduce typos in the numerical code. The 
numerical examples of Section |IV| are excellent tests of 
the code. In Appendix |A| we describe another analytical 
test that we have found very useful to check our imple- 
mentation of curvature quantities. 

As a final comment we note that the condition ^ de- 
termines the time evolution of the determinant 7 of the 
conformally related metric, but not its initial value. The 
latter can be chosen freely in this scheme, in particu- 
lar it does not need to be chosen equal to that of the 
background metric 7 (unlike in the original BSSN formu- 
lation). For some of our numerical simulations, however, 
in particular for the rotating star simulations of Section 
|IV B 2\ we found that rescaling the conformally related 
metric so that its determinant becomes 7 = r 4 sin 2 9 im- 
proved the stability of the simulation, so that it required 
a smaller coefficient 77 in the Kreiss-Oliger dissipation 
term (34) below. 



III. NUMERICAL IMPLEMENTATION 

A. PIRK methods 

The origin of the numerical instabilities in curvilinear 
coordinate systems are related to the presence of stiff 
source terms in the equations, e.g. factors of 1/r 2 oi- 
l/sin 2 (9) that become arbitrary large close to the ori- 
gin or the axis. In the following we will refer to these 
terms as "singular terms". PIRK methods evolve all 
other, i.e. regular, terms in the evolution equations ex- 
plicitly, and then use these updated values to evolve the 
singular terms implicitly. This strategy implies that the 
computational costs of PIRK methods are comparable 
to those of explicit methods. The resulting numerical 
scheme does not need any analytical or numerical inver- 
sion, but is able to provide stable evolutions due to its 
partially implicit component. We refer to [2 6) for a de- 
tailed derivation of PIRK methods (up to third order), 
and limit our discussion here to a simple description of 
the second-order PIRK method that is implemented in 
our code. 

Consider a system of partial differential equations 



u t = d(u,v), 

v t = £ 2 {u) + C 3 {u,v), 



(28) 



where C%, C 2 and £3 are general non-linear differential 
operators. We will denote the corresponding discretized 
operators by Li, L 2 and L 3l respectively. We will further 
assume that L\ and L 3 contain only regular terms, and 
hence will update these terms explicitly, whereas the L 2 
operator contains the singular terms and will therefore 



be treated partially implicitly. Note that L 2 is assumed 
to depend on u only. In the case of the BSSN equations 
this holds for almost all variables; the one exception can 
be treated as discussed in the paragraph below equation 
( B6 ) in Appendix |Bj where we provide the exact form of 
the source terms. 

In our second-order PIRK scheme we update the vari- 
ables it and v from an old timestep n to a new timestep 
n + 1 in two stages. In each of these two stages, we first 
evolve the variable u explicitly, and then evolve the vari- 
able v taking into account the updated values of u for the 
evaluation of the singular L 2 operator. For the system of 
equations ( 28 1 , the first stage 



|W = u n + AtL 1 (u n ,v n ), 



= v" + At 
is followed by the second stage 



^L 2 (u n ) + h 2 (u^)+L 3 (u n ,v n ) 



(29) 



u n+1 = 



,71+1 



1 



+ u {1) +AtL 1 (u^ 1 \v w ) 



At 



n + — [L 2 {u n ) + L 2 {u n+1 ) 

* +L 3 (u n ,v n ) + L 3 (u^ 1 v^) 

(30) 

In the first stage, u is evolved explicitly; the updated 
value it*- 1 ) is used in the evaluation of the L 2 operator for 
the computation of v^K In the second stage, u is again 
evolved explicitly, and the updated value u n+1 is used in 
the evaluation of the L 2 operator for the computation of 
the updated values v n+1 . 

Our PIRK method is stable as long as the timestep is 
limited by a Courant condition; sec cq. ( 33 ) below. 



We include all singular terms appearing in the sources 
of the equations in the L 2 operator. Firstly, the confor- 
mal metric components, hij, the conformal factor, <f>, the 
lapse function, a, and the shift vector, j3 l , are evolved 
explicitly (as u is evolved in the previous PIRK scheme) ; 
secondly, the traceless part of the extrinsic curvature, ay, 
and the trace of the extrinsic curvature, K, are evolved 
partially implicitly, using updated values of a, j3 l , <j> and 
hij] then, the A 1 are evolved partially implicitly, using 
the updated values of a, /3 l , <fi, hij, and K. Finally, 
B l is evolved partially implicitly, using the updated val- 
ues of the previous quantities. Lie derivative terms and 
matter source terms are always included in the explicitly 
treated parts. In Appendix [Bj we give the exact form of 
the source terms included in each operator. 



B. Numerical grid 

We adopt a centered, fourth-order finite differencing 
representation of most spatial derivatives. For each grid 
point, the finite-differencing stencil therefore involves the 
two nearest neighbors in each direction (see Fig. [I]). 
An exception from our fourth-order differencing are ad- 
vective derivatives along the shift, for which we use a 



() 




FIG. 1: A schematic representation of our cell-centered grid 
structure in spherical polar coordinates, for one fixed value 
of (j>. Grid points, marked by the crosses, are placed at the 
center of grid cells, so that no grid point ends up at the center 
(r = 0) or on the axes (sin# = or sin^ = 71"). Our interior 
grid, bordered by solid lines in the figure, covers the region 
< r < r max and < 9 < n (as well as < < 2-71"). As 
suggested by the two highlighted stencils, our fourth-order 
differencing scheme requires two levels of ghost zones outside 
of the interior grid, indicated by the dotted lines. 



second-order (one-sided) upwind scheme. Because of the 
second-order time evolution, and the second-order advec- 
tive terms, our scheme is overall second-order accurate, 
even though for some cases with vanishing shift we have 
found that the error appears to be dominated by the 
fourth-order terms. 

We adopt a cell-centered grid, as shown schematically 
in Fig. [T] Specifically, we divide the physical domain 



covered by our grid, < r < 



< 9 < 7r and < 



cf> < 2tt into N r x N0 x N$ cells with uniform coordinate 
size 

Ar = r max /N r , A6 = 7r/N e , A</> = Tur/N^. (31) 

Because of our fourth-order finite differencing scheme we 
need to pad the interior grid with two layers of ghost 
zones. Except at the outer boundary, each ghost zone 
corresponds to some other zone in the interior of the grid 
(with some other value of 9 and cf>), so that these ghosts 
zones can be filled by copying the corresponding values 
from interior grid points. 

As a concrete example, consider a grid point with an- 
gular coordinates 9 and <f>, say, in the innermost radial 
zone (highlighted by a (blue) filled circle in Fig. [I]) . Eval- 
uating the partial derivative with respect to r at this 
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TABLE I: Parity conditions for components of vectors and 
tensors as implemented in our coordinate-based code. Com- 
ponents of vectors and tensors have to be multiplied with the 
corresponding sign when they are copied into ghost zones at 
the center or the axis. 



point requires two grid points that, formally, have neg- 
ative radii. We can fill these two required ghost points 
by finding the corresponding points in the interior of the 
grid, which have angular coordinates tt — 9 and <f> + w. 
Similarly, evaluating a derivative with respect to 9 for 
a point with angular coordinates (6, (f>) next to the axis 
(see the (red) filled square in Fig. [T]) requires ghost points 
that can be filled by finding the corresponding grid points 
with azimuthal angle <p + ir in the interior of the grid. 

For scalar functions the corresponding function values 
can be copied immediately, but for components of vec- 
tors or tensors, expressed in spherical polar coordinates, a 
possible relative sign has to be taken into account. Essen- 
tially, this occurs because, in spherical polar coordinates, 
the unit vectors may point into the opposite physical di- 
rection when we identify a ghost zone with an interior 
point, i.e. when we go from (9, (j)) to (n — 9, <j) + n) or 
(6, 4> + tt). We list these relative sign changes, as imple- 
mented in our coordinate-based code, in Table [TJ 

We also require two sets of two ghost zones for 0, which 
can be filled directly using periodicity. 

At the outer boundary we also require two ghost zones, 
as suggested by the (red) squared stencil in Fig. [I] We 
impose a Sommerfeld boundary condition, which is an 
approximate implementation of an outgoing wave bound- 
ary condition, to fill these ghost zones. In our coordinate- 
based code we implement this condition by tracing an 
outgoing radial characteristic from each of the outer 
boundary grid points back to the previous time level. 
We then interpolate the corresponding function to the 
intersection of that characteristic and the previous time 
level, and copy that interpolated value, multiplied by a 
suitable fall-off in r, into the boundary grid point. We 
assume a fall-off with r~ 1 for all metric variables (i.e. hij, 
dij, 4> and K) as well as the lapse a, but a r~ 2 fall-off for 
the shift /3 1 as well as X 1 . 

The PIRK method of Section |III A is stable as long 
as the time step At is limited by a Courant-Friedrichs- 
Lewy condition. In order to evaluate this condition we 
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any two grid-points in our cell-centered, spherical polar 
grid. This minimum distance is approximately 



min(Ar, (Ar/2)A0, (Ar/2) sin(A6>/2)A</>). 



We then set 



At = CA n 



(32) 



(33) 



where we have chosen a Courant factor C = 0.4 for all 
simulations in this paper. It is a well-known disadvantage 
of spherical polar coordinates that the accumulation of 
gridpoints in the vicinity of the origin leads to a very 
severe limit on the timestep. We will discuss this issue 
in greater detail in Section [V] 

We use Kreiss-Oliger |30j dissipation to suppress the 
appearance of high frequency noise at late times. Specif- 
ically, we add a term of the form 



iko = 



n 



16At 



(Ar)- 



d 4 / 
dr 4 



(A0) 



(A4>y 



5V 

84> A 



( 34 ) 

to the right hand side of the evolution equation for each 
dynamical variable /. Here r\ is a dimensionless coeffi- 
cient which we have chosen between (for some of our 
short-time evolutions) and 0.001 for the rotating neutron 
star simulation in Sect. IIV B 21 



IV. NUMERICAL EXAMPLES 

A. Weak gravitational waves 

As a first test of our codes we consider small-amplitude 
gravitational waves on a flat Minkowski background. Fol- 
lowing Teukolsky |31j we construct an analytical, linear 
solution for quadrupolar (i — 2) waves from a function 



F(r,t) = A{r±t)e-^ r±t f/ x2 



(35) 



where the constant A is related to the amplitude of the 
wave and A to its wavelength (see also Section 9.1 in 
[14]). We set A = 1, by which all length scales become 
dimensionless. We will consider axisymmetric (to = 0) 
and non-axisymmetric (to = 2) modes separately. 



1. Axisymmetric waves 

We first consider axisymmetric to = waves. Since 
these solutions are independent of the coordinate <p, we 
may choose as small as possible (which is = 2 
in our code) without loss of accuracy. We also choose a 
small amplitude of A — 10~ 7 , so that deviations from the 
analytic solution, which is accurate only to linear order 
in A, are dominated by our finite-difference error, and 
not by terms that are higher-order in A. 

In the following we show results for a numerical grid 
with (AON, 10N, 2) grid points, where N = 1, N = 2 or 
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FIG. 2: Snapshots of the metric coefficient h rr for an axisym- 
metric m = small-amplitude gravitational wave at different 
instances of time. For this simulation we used a grid of size 
(160,40,2) and imposed the outer boundary at r max = 8.0. 
We show data as a function of r in the (arbitrary) direction 
6 — 1.61 and <f> = 4.71. Differences between the numerical 
results (marked by crosses) and the analytical solution (solid 
lines) are smaller than the width of the lines in this graph. 



N = 4, and imposing the outer boundary at r = 8.0. 
For these simulations we used the 1+log lapse condition 
(15), but chose a vanishing shift f3' 1 = instead of the 

161). 



Gamma-driver condition 

In Fig. [2] we show snapshots of the metric function h r 
at different instances of time for our highest-resolution 
simulation with N = 4. For each time, we include the 
numerical results as crosses, as well as the analytical so- 
lution as a solid line. The differences between the nu- 
merical results and analytical solution are well below the 
resolution limit of this graph, so that the two cannot be 
distinguished in this Figure. 

In Fig. [3] we show a convergence test for these waves. 
Specifically, we compute the L 2 -norm of the difference 
between the analytical solution h rr and the analytical 
solution, 



lA/wll = 



(K 



-h 



ana\ 
rr ) 



■dV 



1/2 



(36) 



where V is the coordinate volume of the numerical grid. 
In Fig. [3] we show these norms as a function of time for 
N = 1, N = 2 and N — 4. The norms are rescaled 
with a factor iV 4 ; the convergence of the resulting error 
curves indicates that, at these early times, the error is 
dominated by the fourth-order differencing of the spatial 
derivatives. In spherical polar coordinates, the Courant 
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FIG. 3: The norm of the error in the quantity h rr as a func- 
tion of time for a small-amplitude, axisymmetric gravitational 
wave. We show results for simulations with a grid of size 
(AON, 10N, 2), for N = 1, N = 2 and N = 4, with the outer 
boundary imposed at r = 8.0. At these early times, the error 
appears to be dominated by the fourth-order differencing of 
the spatial derivatives. 



FIG. 4: Snapshots of the metric coefficient h rr for a non- 
axisymmetric m = 2 small-amplitude gravitational wave at 
different instances of time. For this simulation we used a 
grid of size (40, 32, 64) and imposed the outer boundary at 
r mal = 4.0; we show data as a function of r in the direction 
6 — 1.62 and <f) = 3.19. Numerical results are marked by the 
crosses, while the analytical solution is shown as the solid line. 



condition ( 33 ) limits the time step to such small values 



that the second-order errors associated with our PIRK 
method are smaller than the fourth-order error of our 
spatial derivatives (for vanishing shift). 



2. Nonaxisymmetric waves 

Non-axisymmetric gravitational waves represent a 
rare example of an analytical, time-dependent, three- 
dimensional, albeit weak-field solution to the Einstein 
equations. Clearly, this solution represents a stringent 
test for our code. 3 

In Fig. g] we show results for an m = 2 wave, again 
for an amplitude A = 10~ 7 . As in Fig. [2J we graph 
solutions for h rr as functions of r at different instances of 
time. Again, our numerical solution (marked by crosses) 
can hardly be distinguished from the analytical solution 
(shown as solid lines). 



B. Hydro-without-hydro 

As a test of strong-field, but regular solutions we con- 
sider spacetimes containing relativistic stars. In general, 
this requires evolving the stellar matter self-consistcntly 
with the gravitational fields, for example by solving the 
equations of relativistic hydrodynamics. Since this is be- 
yond the scope of this paper, we here adopt the "hydro- 
without-hydro" approach suggested by [35]. In this ap- 
proach, which can also be described as an "inverse- 
Cowling approximation" , we leave the matter sources 
fixed, and evolve only the gravitational fields. In this 
way, it is possible to assess the stability of a spacetime 
evolution code, and its capability of accurately evolving 
strong but regular gravitational fields in spacetimes with 
static matter, without having to worry about the hydro- 
dynamical evolution. These simulations serve as both a 
testbed and a preliminary step towards fully relativistic 
hydrodynamical simulations of stars. In this Section we 
consider static and uniformaly rotating stars separately. 



For a code in Cartesian coordinates, even a spherically symmetric 
solution represents a stringent test, because the symmetry is not 
reflected by the numerical grid. In our code, however, numerical 
expressions simplify for spherical or axisymmetric solutions, so 
that they do not test every aspect of the code. 



1. Spherical neutron stars 

We first consider non-rotating relativistic stars, de- 
scribed by the Tolman-Oppenheimer-Volkoff (TOV) so- 
lution [331 [31] . We focus on a polytropic TOV star with 
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FIG. 5: Snapshots of the conformal exponent cj> and the lapse 
a for a rapidly rotating star (see text for details). We show 
both functions both at the initial time, and at a late time 
t — 318M. We also show both functions along rays in two 
different directions, one very close to the equator, the other 
pointing close to the pole. Both profiles remain very similar 
to their initial data throughout the evolution. 



polytropic index T = 2, and with a gravitational mass 
of about 85% of the maximum-allowed mass. For this 
model, the central density is about 40% of that of the 
maximum mass model. We evolved this star with the 
1+log slicing condition for the lapse (15), but kept the 
shift fixed to zero. Because the spacetime is spherically 
symmetric, we could choose both Ng and N$ as small as 
possible (Ng = = 2) without loss of accuracy. 



Even for very modest grid resolutions in the radial di- 
rection (e.g. N r = 40, with the outer boundary imposed 
at four times the stellar radius), we found that the grav- 
itational fields settle down into an equilibrium that is 
similar to the initial data. After this initial transition, 
which is caused by the finite-difference error, the stel- 
lar surface as well as the outer boundaries (see [12])) the 
solution remains stable. 



2. Rotating neutron stars 

The evolution of the spacetime of a rapidly rotating rel- 
ativistic star is a more demanding test than the previous 
one, as it breaks spherical symmetry and instead involves 
axisymmetric non- vacuum initial data in the strong grav- 
ity regime. The initial data used for this test are the nu- 
merical solution of a stationary and axisymmetric equi- 
librium model of a rapidly and uniformly rotating rel- 



ativistic star [35], which is computed using the Lorene 
code [36] . 

We consider a uniformly rotating star with the same 
r = 2 polytropic equation of state as the non-rotating 
model of Sect. |IVB~T Our particular model has the same 
central rest-mass density as that non-rotating model, but 
rotates at 92% of the allowed mass-shedding limit (for a 
star of that central density); expressed in terms of the 
gravitational mass M, the corresponding spin period is 
approximately 157 M. The ratio of the polar to equa- 
torial coordinate radii for this model is 0.7. For this 
simulation we adopted both the 1+log condition for the 
lapse ( 15 ) and the Gamma-driver condition for the shift 

For this test we adopted a grid of size (48,32,2), and 
imposed the outer boundary at 25. 5M, which equals four 
times the equatorial radius. In Fig. [5] we show the initial 
and late-time profiles of the conformal exponent <j) and 
the lapse a, both in a direction close to the equator and 
close to the axis. Evidently, both functions remain very 
close to their initial values throughout the evolution, as 
they should. 



C. Schwarzschild 

In this Section we present results for two different sim- 
ulations involving Schwarzschild black holes. In Section 
IIVC II we evolve a Schwarzschild black hole in a "trum- 
pet" geometry [37U39] . which, in the limit of infinite res- 
olution, is a time-independent solution to the Einstein 
equations given our slicing conditions (15). In Section 



IV C 2 we adopt wormholc initial data and follow the co- 



ordinate transition to a trumpet geometry. 



1. Trumpet initial data 

Maximally sliced trumpet data [38] represent a time- 
independent slicing of the Schwarzschild spacetime that 
satisfies our slicing condition (15). The solution can 



be expressed analytically in isotropic coordinates, albeit 
only in parametrized form |40] . In this Section we adopt 
these trumpet data as initial data, so that, in the con- 
tinuum limit, the solution should remain independent of 
time. 

For trumpet data the conformal factor diverges at 
r = 0. While, on our cell-centered grid, functions are 
never evaluated directly at the origin, derivatives in the 
neighborhood of the singularity at the origin are clearly 
affected by the singular behavior of the conformal fac- 
tor. However, the great virtue of the "moving-puncture" 
gauge conditions |l5| and ( fl6] ) is that these errors only 
affect the neighborhood of the puncture, and do not spoil 
the evolution globally (2] [3] |37l [41] . In the following we 
will demonstrate these properties in our code using spher- 
ical polar coordinates. 
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FIG. 6: The difference between the maximum of the radial 
component of the shift vector, /3maxi an d its initial value 
/3max|t=o, as a function of time. For these simulations we 
used a grid of size (160iV,2,2) for N = 1, 2, 4 and 8, and 
imposed the outer boundary at r max = 16. 0M. We rescale 
all differences with iV 2 , so that the convergence of these lines 
demonstrates second-order convergence. 



For the simulations presented in this section we 
adopted a numerical grid of size size (1607V, 2, 2) for 
N = 1, 2, 4 and 8, with the outer boundary imposed 
at r max = 16.0A/ . 

In Fig. [6] we show results for the maximum of the radial 
component j3 r of the shift vector as a function of time. 
Specifically, we show the difference between these maxi- 
mum values /3^ ax and their initial values /3maxl«=o- Since 
our trumpet data represent a time-independent solution 
to the Einstein equations and our slicing and gauge con- 
ditions (15) and (16), these differences should converge 
to zero as the grid resolution is increased. In Fig. [6] we 
multiply the differences with iV 2 ; the convergence of the 
resulting lines therefore demonstrates second-order con- 
vergence of the simulation. Apparently the error in these 
simulations is dominated by the second-order advective 
terms. 



We also note that the outer boundary introduces error 
terms that depend on both the grid resolution and the 
location of the outer boundary. Since the latter does not 
decrease when we increase the grid resolution, the code 
converges more slowly in regions that have come into 
causal contact with the outer boundary. We therefore 
include in Fig. [6] only sufficiently early times, before the 
location of the shift's maximum is affected by the outer 
boundary. 



FIG. 7: Initial data and final profiles of the conformal ex- 
ponent <f>, the lapse function a, and the shift /3 r , showing 
the coordinate transition from wormhole initial data to time- 
independent trumpet data. The (blue) long-dashed lines rep- 
resent the initial data at t = 0, the (red) dashed lines show our 
numerical results at time t = 79M, and the (black) solid lines 
show the analytical trumpet solution [40]. The initial data 
appear double-valued because we graph this functions as a 
function of the areal radius R (see text for details). For these 
simulations we adopted a grid size (10240, 2, 2) with the outer 
boundary imposed at r = 256M. (In these graphs we did not 
include the innermost two grid points, which are affected by 
the singular behavior of the puncture.) 



2. Wormhole initial data 

We now turn to evolutions of wormhole initial data, 
representing a horizontal slice through the Penrose dia- 
gram of a Schwarzschild black hole. For these data, the 
conformal factor is given by 



ip = 1 + 



M 
27' 



(37) 



the conformally related metric is flat, 7^ = rjij, and the 
extrinsic curvature vanishes, Aij = = K. Instead of 
choosing the Killing lapse and Killing shift, which would 
leave these data time-independent, we choose, at t = 0, 
a "pre-collapsed" lapse [8] 



a = ip 



(38) 



and a vanishing shift, f3 l = 0. We then evolve the 



lapse and the shift with the 1+log condition ( 15 ) and 
the Gamma-driver condition (16). 

Since these initial data do not represent a time- 
independent solution to the Einstein equations together 
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FIG. 8: The maximum of the radial shift /3 r as a function 
of time. We show results for different grid sizes (1280TV, 2, 2) 
for TV = 1, 2, 4 and 8, with the outer boundary imposed at 
256M. After a brief transition from the initial data /?'" = 0, 
the shift settles down into a new equilibrium. For relatively 
coarse grid resolutions the shift experiences a slow drift, but 
this drift disappears as the grid resolution is increased. 



FIG. 9: Profiles of violations of the Hamiltonian constraint 
Q at time t = 79M. As in Fig. [8] we show results for 
grid sizes (1280TV, 2, 2) for TV = 1, 2, 4 and 8, with the outer 
boundary imposed at 256M. All results are rescaled with TV 2 ; 
the convergence of the resulting lines demonstrates second- 
order convergence of our code. 



with our gauge conditions, we observe a non-trivial time 
evolution that represents a coordinate evolution. For 
the "non-advective" 1+log condition ( [15] ), this coordi- 
nate transition results in the maximally sliced trumpet 
geometry of Section IV C 2 39 . In Fig. [7] we show this 
coordinate transition for the conformal exponent <f>, the 
lapse a and the shift f3 r . 

We note that some care has to be taken when the nu- 
merical and analytical results are compared. The an- 
alytical solution of [30] assumes 7-y = r)ij. We also 
choose 7y = rjij in our initial data, but this relation 
is not necessarily maintained during the time evolution, 
so that the numerical and analytical solutions may be 
represented in different spatial coordinate systems (but 
on the same spatial slice). In order to compare the two 
solutions we therefore graph all quantities as a function 
of the gauge-invariant areal radius R. Since for wormhole 
data each value of R > 2M corresponds to two values of 
the isotropic radius r, the initial data in Fig. [7] appear 
double- valued. For these comparisons with the analyt- 
ical solution we also graph the orthonormal component 
of the shift /3 r rather than the coordinate component f3 r 
itself. Fig. [7] clearly shows the coordinate transition from 
wormhole initial data to the trumpet equilibrium solu- 
tion. 

In Fig. [8] we show the maximum of the radial shift /3 r 
as a function of time. After a brief period of a coordinate 



transition the shift settles down into a new equilibrium. 
We show results for grid sizes (1280TV, 2, 2) for TV = 1, 2, 
4 and 8, with the outer boundary imposed at 256M. The 
graph shows that differences between the different results 
decrease rapidly as the grid resolution is increased. For 
our coarser grid resolutions the shift still experiences a 
slow drift after the initial transition, but this drift de- 
creases as the grid resolution is increased. 

Finally, in Fig. [9] we show profiles of the violations of 
the Hamiltonian constraint ( 13 ) at time t — 79M. In this 
graph all results are rescaled with TV 2 ; the convergence of 
the resulting lines demonstrates that the numerical error 
in these simulations is again dominated by the second- 
order implementations of the advective shift terms, and 
possibly the time evolution. 



V. DISCUSSION 

In this paper we demonstrate that a PIRK method can 
be used to solve the Einstein equations in spherical po- 
lar coordinates without any need for any regularization 
at the origin or on the axis. Specifically, we integrate 
a covariant version of the BSSN equations in three spa- 
tial dimensions without any symmetry assumptions. To 
the best of our knowledge, these calculations represent 
the first successful three-dimensional numerical relativity 
simulations using spherical polar coordinates. We con- 
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sider several test cases to assess the stability, accuracy 
and convergence of the code, namely weak-held "Teukol- 
sky" gravitational waves, "hydro-without-hydro" simula- 
tions of static and rotating relativistic stars, and single 
black holes. 

Spherical polar coordinates have several advantages 
and disadvantages over Cartesian coordinates. At least 
in single-grid calculations, spherical polar coordinates al- 
low for a more effective allocation of the numerical grid 
points for applications that involve one center of mass, 
for example gravitational collapse of single stars or su- 
pernovae. This is true even for uniform grids, which we 
adopt in this paper, but curvilinear coordinate systems 
also facilitate the use of non-uniform grids (e.g. a log- 
arithmic radial coordinate) to achieve a high resolution 
near the origin while keeping the outer boundary suffi- 
ciently far. 

Spherical polar coordinates have another strong advan- 
tage over Cartesian coordinates. In simulations of super- 
novae or gravitational collapse, for example, the shape of 
the stellar objects is not well represented by Cartesian 
grids. This mismatch between the symmetry of the ob- 
ject and the grid creates direction-dependent numerical 
errors, which are observed to trigger m = 4 modes that 
grow in time. Since spherical polar coordinates mimic 
the symmetry of collapsing stars more accurately, we ex- 
pect that this problem can at least be reduced with these 
coordinates. 

However, spherical polar coordinates also have disad- 
vantages. One of these disadvantages is of practical na- 
ture: the equations in spherical polar coordinates include 
many more terms than those in Cartesian coordinates, 
which makes the numerical implementation more cum- 
bersome and error prone. Spherical polar coordinates 
also introduce coordinate singularities that traditionally 
have created many numerical problems; but these prob- 
lems can be avoided when using a PIRK method. 

Perhaps the most severe disadvantage of spherical po- 
lar coordinates is caused by the Courant limitation on the 
time step. As shown in eq. (33), the close proximity of 



pletely implicit scheme, so that the time step there is no 
longer limited by the Courant condition ( 33 1 . Similar 



grid points close to the origin limits the size of the time 
steps At to increasingly small values as the resolution 
is increased. In three-dimensional simulations, At de- 
creases approximately with the product At oc ArAQAcf). 
This is a severe disadvantage compared to Cartesian co- 
ordinates where typically At oc Ax 1 . However, this prob- 
lem is not unique to numerical relativity, and instead is 
well-known from dynamical simulations in spherical po- 
lar coordinates in any held. Accordingly, several differ- 
ent approaches to either solving or reducing this problem 
have been suggested. 

One possible approach is to reduce the grid resolution 
in the angular directions, Ng and JV^, close to the ori- 
gin. However, for many applications the angular depen- 
dence of the solution may be independent of the radius, 
so that this approach might severely limit the accuracy 
of the results. It may also be possible to replace the 
PIRK method in a sphere around the origin with a corn- 



implicit /explicit (IMEX) "split-by-region" schemes have 
been suggested, for example, in [35] in the context of 
spectral schemes. Finally, the "Yin- Yang" method sug- 
gested in [43] mitigates the restrictions imposed by the 
Courant condition (33 1 as follows. Note that the small- 



est physical distance between grid points, which in turn 
limits the time step At, occurs next to the axis. In the 
Yin- Yang method, the unit sphere is therefore covered 
by two different grids that are rotated by an angle of 90 
degrees with respect to each other. Each one covers only 
a region around its equator, thereby avoiding the most 
severe limitation on the time step next to the axis, but 
combined both grids cover the entire unit sphere. 

Despite the small time step, however, we have been 
able to complete all simulations presented in this paper 
even with a serial code - in fact, some of our simulations 
were performed on a laptop computer. 
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Appendix A: A numerical test for curvature 
quantities in spherical polar coordinates 

In spherical polar coordinates, in particular in the ab- 
sence of any symmetry assumptions, the numerical im- 
plementation of curvature quantities involves a signifi- 
cant number of terms that can easily introduce mistakes 
(see Section II B ) . One way of testing this part of the nu- 
merical code is to compare with known analytical solu- 
tions, for example for the Schwarzschild metric. However, 
most analytical solutions feature symmetries (e.g. spher- 
ical symmetry for Schwarzschild) that simplify the prob- 
lem in the spherical polar coordinates of our code. As 
a consequence, many terms vanish identically for these 
solutions, so that not all terms in the code are tested. In 
this Appendix we describe a simple test that is also ana- 
lytical, but is neither spherically nor axially symmetric, 
and hence a very stringent test. 

Starting with the hat metric in Cartesian coordinates 
we introduce a coordinate transformation of each coordi- 
nate x l that only depends on that coordinate itself; the 
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resulting metric then takes the form 

f{x) 
= I g(y) 
h{z) 



(Al) 



where fix), giy) and /i(z) are arbitrary functions. Trans- 
forming this metric into spherical polar coordinates leads 
to a metric for which, in general, all coefficients are non- 
zero and depend on the coordinates in potentially com- 
plicated ways. 

In Cartesian coordinates, the only non-vanishing 
Christoffel symbols are 
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where the prime denotes a derivative with respect to the 
argument. Given that the Ar* fc transform like tensors, 
we can obtain the corresponding coefficients in spherical 
polar coordinates with a simple coordinate transforma- 
tion. For sufficiently general functions f(x), giy) and 
h(z), all 18 components of APy in spherical polar coor- 
dinates will be non-zero. This yields analytical expres- 
sions for the connection coefficients (23 1 that can then 
be compared with numerical results. 



Similarly, the connection functions are given by 

'fix) g'(y) hf^y 
2/ 2 ' 2g 2 ' 2h 2 



3.4mA l = 



(A5) 



in Cartesian coordinates, and can be transformed into 
spherical polar coordinates with a simple coordinate 
transformation. 

Finally, all components of the Ricci tensor in spheri- 
cal polar coordinates should converge to zero, since the 
metric (All is still flat. In Fig. 10 we show numerical 



examples for 



m 

g{y) 

Kz) 



0.1 x 2 
0.3 y 2 
0.5 z 2 
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All components of Rij are non-zero, but converge to zero 
as the grid resolution is increased. In the graph we rescale 
all results with TV 4 , so that the convergence of the result- 
ing quantities indicates fourth-order convergence of our 
implementation of the Ricci tensor, as expected. 



Appendix B: Detailed source terms included in the 
PIRK operators for the evolution equations 
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FIG. 10: Values of the norms of th e co mponents of the Ricc i 
tensor, H-RijH, for the flat metric (All with functions (A6l, 
evaluated using grid sizes (16 A, 8 A, 167V) for A = 1, 2, 4 and 
8. All values are rescaled with A 4 , so that the convergence of 
these results indicates fourth-order convergence of our imple- 
mentation. 



equations into the explicit and partially implicit opera- 
tors. 

We start each time step by evolving the conforma! met- 
ric components, hy, the conformal factor 0, the lapse 
function, a, and the shift vector, /3 J , explicitly, i.e., all 
the source terms of the evolution equations of these vari- 
ables are included in the L\ operator of the second-order 
PIRK method. 

We then evolve the traceless part of the extrinsic cur- 
vature, CLij, and the trace of the extrinsic curvature, K, 
partially implicitly. More specifically, the corresponding 
Li and L3 operators associated with the evolution equa- 
tions for a>ij and K in terms of t he o riginal BSSN variable 
Aij, related to 



through cq. (21), are 
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+ AD^aD 



3Y 
TF 



= --AijD k j3 k - 2aA ik A k j + aA tJ K, 



J 3(A zj ) - 3 ~ n 



(Bl) 
(B2) 



2(K) 
3(K) 
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iD 2 a + 2D i aD i $) + a Aij A 13 , (B3) 

(B4) 



9}, ([T5|-([l6| using a 



We evolve the evolution eqs 
second-order PIRK method. In this Appendix we pro- 
vide details on how we split the right-hand sides of these 



The A 1 are evolved partially implicitly, using the updated 
values of a, /3 4 , 4>, hij, and K. In terms of the original 
BSSN variable A\ related to A* through eq. |22l, the 



14 



operators are 



L 



2(A< 



l . -r —D l Djfij 



2A jk (5 i j d k a-6a5 l j d k (j) 
2 

: —i 

3 



-aj ij djK 



(B5) 
(B6) 



We note that the evaluation of the Ricci tensor i?,-,- in 



availab 



eq. flBl I requires updated values of A* before they become 



It is possible to either replace these updated 



values with old values, or to update the A 1 provisionally 



in a purely explicit step, use these values in eq. (Bl|, 



but then overwrite these values after the A 1 are updated 
partially implicitly. We have used the latter approach in 
the simulations presented in this paper. 

Finally, the B l are evolved partially implicitly, using 
the updated values of the previous quantities, 



L 



2(s*) - ^A 4 , 



(B7) 
(B8) 



Matter source terms and Lie derivative terms are always 
included in the explicitly treated parts. 
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